Purpose:

This is the final workflow file. We take processed data from the FinalData folder and do some point pattern analysis with the pokemon point level data.

It seems that pokemon point spawn locations are not so much determined by neighborhood characteristics, but by where other points are! Let’s investigate the tendency of pokemon spawn points to co-locate further, through point pattern analysis (distance to nth nearest neighbors).

0 Libraries

library(spatstat)
library(sf)
library(tidyverse)
setwd("~/Desktop/GIS3FinalProject")

Step 1: Pre processing

1.1 Load point data in sf (epsg 7131)

poke_points  <- 
  st_read("./FinalData/ShapefileVersions/NoNAPokePointsFinal7131.shp")  
Reading layer `NoNAPokePointsFinal7131' from data source `/Users/calvinzhang/Desktop/GIS3FinalProject/FinalData/ShapefileVersions/NoNAPokePointsFinal7131.shp' using driver `ESRI Shapefile'
Simple feature collection with 9527 features and 21 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 41393.58 ymin: 19369.43 xmax: 59579.32 ymax: 35780.02
Projected CRS: Transverse_Mercator
plot(poke_points)
plotting the first 9 out of 21 attributes; use max.plot = 21 to plot all

1.2 Convert points shapefile to spatstat format (ppp)

poke_points  <- as.ppp(poke_points)
only first attribute column is used for marks
marks(poke_points) <- NULL
poke_points.km <- rescale(poke_points, 1000, "km") # EPSG 7131 is in meters, times 1000 to km
plot(poke_points.km)

Step 2: Analyze Average distance to the nth nearest neighbor point

2.1: Calculate average distance to the nth nearest neighbor point

mean(nndist(poke_points.km, k=1))
[1] 0.05224023
mean(nndist(poke_points.km, k=2))
[1] 0.08158853
# ANN stands for average nth neighbor (ANN) distance
ANN <- apply(nndist(poke_points.km, k=1:500),2,FUN=mean)

2.2: Plot average distance to the 1st, 2nd, …, 500th nearest neighbors

n_nearest_dist <- data.frame(nth_nearest = 1:500, distkm = ANN)

ggplot(data = n_nearest_dist) +
  geom_point(aes(x=nth_nearest, y=distkm,col=-distkm)) +
  ggtitle("Distance in km to nth nearest neighbor") +
  xlab("Nth nearest neighbor") + 
  ylab("Average distance to neighbor in km") + 
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5))

Step 3 Comparing actual 1st nearest neighbor distances to randomized results

3.1 Create randomnized points and calculate their average 1st nearest neighbor distances

# Simulate random points N times, and for each simulation, store the average 1st nearest neighbor distances in ANN.r

n     <- 599L               # Number of simulations
ANN.r <- vector(length = n) # Vector to store random (r) average nearest neighbor distances (ANN)
for (i in 1:n){
  rand.p   <- rpoint(n=poke_points.km$n, win=poke_points.km)  # Generate random point locations
  ANN.r[i] <- mean(nndist(rand.p, k=1))  # Tally the average neearest neighbor distance
}

plot(rand.p)

3.2 Compare randomnized points’ average nearest neighbor distances to actual average nearest neighbor distance

# Actual distance
ANN.p <- mean(nndist(poke_points.km, k=1))

# Convert vector of simulated average nearest neighbor distances to dataframe
random_nearest1 <- data.frame(random_dist = ANN.r)

# Plot results
ggplot(data = random_nearest1) +
  geom_histogram(aes(x=random_dist), fill="blue", bins=2000) +
  geom_vline(aes(xintercept = ANN.p), col='red') +
  xlab("Distance to 1st nearest neighbor, km") +
  ylab("Frequency") + 
  ggtitle("Frequency of distance in km to nth nearest neighbor, 
          actual (red) vs.random (blue) ") +
  theme(legend.position = "none", plot.title = element_text(hjust = 0.5)) 

3.3 Calculate pseudo p value for significance of average 1st nearest neighbor distance

N.greater <- sum(ANN.r > ANN.p)
pseudo_p <- min(N.greater + 1, n + 1 - N.greater) / (n +1)
pseudo_p
[1] 0.001666667
LS0tCnRpdGxlOiAiTmVhcmVzdCBuZWlnaGJvciBhbmFseXNpcyBvZiBQb2tlbW9uIHBvaW50cyIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQojIFB1cnBvc2U6IApUaGlzIGlzIHRoZSBmaW5hbCB3b3JrZmxvdyBmaWxlLiBXZSB0YWtlIHByb2Nlc3NlZCBkYXRhIGZyb20gdGhlIEZpbmFsRGF0YSBmb2xkZXIgYW5kIGRvIHNvbWUgcG9pbnQgcGF0dGVybiBhbmFseXNpcyB3aXRoIHRoZSBwb2tlbW9uIHBvaW50IGxldmVsIGRhdGEuCgpJdCBzZWVtcyB0aGF0IHBva2Vtb24gcG9pbnQgc3Bhd24gbG9jYXRpb25zIGFyZSBub3Qgc28gbXVjaCBkZXRlcm1pbmVkIGJ5IG5laWdoYm9yaG9vZCBjaGFyYWN0ZXJpc3RpY3MsIGJ1dCBieSB3aGVyZSBvdGhlciBwb2ludHMgYXJlISBMZXQncyBpbnZlc3RpZ2F0ZSB0aGUgdGVuZGVuY3kgb2YgcG9rZW1vbiBzcGF3biBwb2ludHMgdG8gY28tbG9jYXRlIGZ1cnRoZXIsIHRocm91Z2ggcG9pbnQgcGF0dGVybiBhbmFseXNpcyAoZGlzdGFuY2UgdG8gbnRoIG5lYXJlc3QgbmVpZ2hib3JzKS4gCgojIDAgTGlicmFyaWVzCmBgYHtyfQpsaWJyYXJ5KHNwYXRzdGF0KQpsaWJyYXJ5KHNmKQpsaWJyYXJ5KHRpZHl2ZXJzZSkKc2V0d2QoIn4vRGVza3RvcC9HSVMzRmluYWxQcm9qZWN0IikKYGBgCgojIFN0ZXAgMTogUHJlIHByb2Nlc3NpbmcKIyMgMS4xIExvYWQgcG9pbnQgZGF0YSBpbiBzZiAoZXBzZyA3MTMxKQpgYGB7cn0KcG9rZV9wb2ludHMgIDwtIAogIHN0X3JlYWQoIi4vRmluYWxEYXRhL1NoYXBlZmlsZVZlcnNpb25zL05vTkFQb2tlUG9pbnRzRmluYWw3MTMxLnNocCIpICAKCnBsb3QocG9rZV9wb2ludHMpCmBgYAojIyAxLjIgQ29udmVydCBwb2ludHMgc2hhcGVmaWxlIHRvIHNwYXRzdGF0IGZvcm1hdCAocHBwKQpgYGB7cn0KcG9rZV9wb2ludHMgIDwtIGFzLnBwcChwb2tlX3BvaW50cykKbWFya3MocG9rZV9wb2ludHMpIDwtIE5VTEwKcG9rZV9wb2ludHMua20gPC0gcmVzY2FsZShwb2tlX3BvaW50cywgMTAwMCwgImttIikgIyBFUFNHIDcxMzEgaXMgaW4gbWV0ZXJzLCB0aW1lcyAxMDAwIHRvIGttCnBsb3QocG9rZV9wb2ludHMua20pCmBgYAojIFN0ZXAgMjogQW5hbHl6ZSBBdmVyYWdlIGRpc3RhbmNlIHRvIHRoZSBudGggbmVhcmVzdCBuZWlnaGJvciBwb2ludAojIyAyLjE6IENhbGN1bGF0ZSBhdmVyYWdlIGRpc3RhbmNlIHRvIHRoZSBudGggbmVhcmVzdCBuZWlnaGJvciBwb2ludApgYGB7cn0KbWVhbihubmRpc3QocG9rZV9wb2ludHMua20sIGs9MSkpCm1lYW4obm5kaXN0KHBva2VfcG9pbnRzLmttLCBrPTIpKQoKIyBBTk4gc3RhbmRzIGZvciBhdmVyYWdlIG50aCBuZWlnaGJvciAoQU5OKSBkaXN0YW5jZQpBTk4gPC0gYXBwbHkobm5kaXN0KHBva2VfcG9pbnRzLmttLCBrPTE6NTAwKSwyLEZVTj1tZWFuKQpgYGAKCiMjIDIuMjogUGxvdCBhdmVyYWdlIGRpc3RhbmNlIHRvIHRoZSAxc3QsIDJuZCwgLi4uLCA1MDB0aCBuZWFyZXN0IG5laWdoYm9ycwpgYGB7cn0Kbl9uZWFyZXN0X2Rpc3QgPC0gZGF0YS5mcmFtZShudGhfbmVhcmVzdCA9IDE6NTAwLCBkaXN0a20gPSBBTk4pCgpnZ3Bsb3QoZGF0YSA9IG5fbmVhcmVzdF9kaXN0KSArCiAgZ2VvbV9wb2ludChhZXMoeD1udGhfbmVhcmVzdCwgeT1kaXN0a20sY29sPS1kaXN0a20pKSArCiAgZ2d0aXRsZSgiRGlzdGFuY2UgaW4ga20gdG8gbnRoIG5lYXJlc3QgbmVpZ2hib3IiKSArCiAgeGxhYigiTnRoIG5lYXJlc3QgbmVpZ2hib3IiKSArIAogIHlsYWIoIkF2ZXJhZ2UgZGlzdGFuY2UgdG8gbmVpZ2hib3IgaW4ga20iKSArIAogIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICJub25lIiwgcGxvdC50aXRsZSA9IGVsZW1lbnRfdGV4dChoanVzdCA9IDAuNSkpCmBgYAojIFN0ZXAgMyBDb21wYXJpbmcgYWN0dWFsIDFzdCBuZWFyZXN0IG5laWdoYm9yIGRpc3RhbmNlcyB0byByYW5kb21pemVkIHJlc3VsdHMKCiMjIDMuMSBDcmVhdGUgcmFuZG9tbml6ZWQgcG9pbnRzIGFuZCBjYWxjdWxhdGUgdGhlaXIgYXZlcmFnZSAxc3QgbmVhcmVzdCBuZWlnaGJvciBkaXN0YW5jZXMKYGBge3J9CiMgU2ltdWxhdGUgcmFuZG9tIHBvaW50cyBOIHRpbWVzLCBhbmQgZm9yIGVhY2ggc2ltdWxhdGlvbiwgc3RvcmUgdGhlIGF2ZXJhZ2UgMXN0IG5lYXJlc3QgbmVpZ2hib3IgZGlzdGFuY2VzIGluIEFOTi5yCgpuICAgICA8LSA1OTlMICAgICAgICAgICAgICAgIyBOdW1iZXIgb2Ygc2ltdWxhdGlvbnMKQU5OLnIgPC0gdmVjdG9yKGxlbmd0aCA9IG4pICMgVmVjdG9yIHRvIHN0b3JlIHJhbmRvbSAocikgYXZlcmFnZSBuZWFyZXN0IG5laWdoYm9yIGRpc3RhbmNlcyAoQU5OKQpmb3IgKGkgaW4gMTpuKXsKICByYW5kLnAgICA8LSBycG9pbnQobj1wb2tlX3BvaW50cy5rbSRuLCB3aW49cG9rZV9wb2ludHMua20pICAjIEdlbmVyYXRlIHJhbmRvbSBwb2ludCBsb2NhdGlvbnMKICBBTk4ucltpXSA8LSBtZWFuKG5uZGlzdChyYW5kLnAsIGs9MSkpICAjIFRhbGx5IHRoZSBhdmVyYWdlIG5lZWFyZXN0IG5laWdoYm9yIGRpc3RhbmNlCn0KCnBsb3QocmFuZC5wKQpgYGAKCiMjIDMuMiBDb21wYXJlIHJhbmRvbW5pemVkIHBvaW50cycgYXZlcmFnZSBuZWFyZXN0IG5laWdoYm9yIGRpc3RhbmNlcyB0byBhY3R1YWwgYXZlcmFnZSBuZWFyZXN0IG5laWdoYm9yIGRpc3RhbmNlCmBgYHtyfQojIEFjdHVhbCBkaXN0YW5jZQpBTk4ucCA8LSBtZWFuKG5uZGlzdChwb2tlX3BvaW50cy5rbSwgaz0xKSkKCiMgQ29udmVydCB2ZWN0b3Igb2Ygc2ltdWxhdGVkIGF2ZXJhZ2UgbmVhcmVzdCBuZWlnaGJvciBkaXN0YW5jZXMgdG8gZGF0YWZyYW1lCnJhbmRvbV9uZWFyZXN0MSA8LSBkYXRhLmZyYW1lKHJhbmRvbV9kaXN0ID0gQU5OLnIpCgojIFBsb3QgcmVzdWx0cwpnZ3Bsb3QoZGF0YSA9IHJhbmRvbV9uZWFyZXN0MSkgKwogIGdlb21faGlzdG9ncmFtKGFlcyh4PXJhbmRvbV9kaXN0KSwgZmlsbD0iYmx1ZSIsIGJpbnM9MjAwMCkgKwogIGdlb21fdmxpbmUoYWVzKHhpbnRlcmNlcHQgPSBBTk4ucCksIGNvbD0ncmVkJykgKwogIHhsYWIoIkRpc3RhbmNlIHRvIDFzdCBuZWFyZXN0IG5laWdoYm9yLCBrbSIpICsKICB5bGFiKCJGcmVxdWVuY3kiKSArIAogIGdndGl0bGUoIkZyZXF1ZW5jeSBvZiBkaXN0YW5jZSBpbiBrbSB0byBudGggbmVhcmVzdCBuZWlnaGJvciwgCiAgICAgICAgICBhY3R1YWwgKHJlZCkgdnMucmFuZG9tIChibHVlKSAiKSArCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gIm5vbmUiLCBwbG90LnRpdGxlID0gZWxlbWVudF90ZXh0KGhqdXN0ID0gMC41KSkgCmBgYAojIyAzLjMgQ2FsY3VsYXRlIHBzZXVkbyBwIHZhbHVlIGZvciBzaWduaWZpY2FuY2Ugb2YgYXZlcmFnZSAxc3QgbmVhcmVzdCBuZWlnaGJvciBkaXN0YW5jZQpgYGB7cn0KTi5ncmVhdGVyIDwtIHN1bShBTk4uciA+IEFOTi5wKQpwc2V1ZG9fcCA8LSBtaW4oTi5ncmVhdGVyICsgMSwgbiArIDEgLSBOLmdyZWF0ZXIpIC8gKG4gKzEpCnBzZXVkb19wCmBgYAoK